Los datos diabetes-data.csv contienen la información de medición de glucosa en la sangre de 70 pacientes con diabetes observados en diferentes momentos del tiempo (horas, días, semanas, meses).
El contenido y formato de los datos es el siguiente:
Los datos de concentración en la sangre se obtienen de dos fuentes:
Las mediciones eletrónicas automáticas se obtienen con etiquetas de tiempo específicas, mientras que las mediciones obtenidas de los registros escritos corresponden a horarios ficticios con la siguiente relación: Desayuno (08:00), comida (12:00), cena (18:00), y hora de dormir (22:00).
El campo code corresponde al tipo de insulina administrada al paciente.
Nota: La información de la adminitración de insulina, en este momento, es irrelevante.
date time code value
Length:29264 Length:29264 Min. : 4.0 Min. : 0.00
Class :character Class1:hms 1st Qu.:33.0 1st Qu.: 6.00
Mode :character Class2:difftime Median :48.0 Median : 24.00
Mode :numeric Mean :46.5 Mean : 79.42
3rd Qu.:60.0 3rd Qu.:142.00
Max. :72.0 Max. :501.00
NA's :8
individual atributo
Min. : 1.00 Length:29264
1st Qu.:21.00 Class :character
Median :34.00 Mode :character
Mean :36.45
3rd Qu.:55.00
Max. :70.00
Observations: 29,264
Variables: 6
$ date <chr> "04-21-1991", "04-21-1991", "04-21-1991", "04-21-19...
$ time <time> 09:09:00, 09:09:00, 09:09:00, 17:08:00, 17:08:00, ...
$ code <int> 58, 33, 34, 62, 33, 48, 58, 33, 34, 33, 62, 33, 58,...
$ value <dbl> 100, 9, 13, 119, 7, 123, 216, 10, 13, 2, 211, 7, 25...
$ individual <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
$ atributo <chr> "medición de glucosa en sangre antes del desayuno",...
A partir de lo anterior, puede extenderse el modelo para considerar un modelo jerárquico que responda si existe o no diferencias entre las mediciones.
El modelo se define de la siguiente manera:
\[ log(y_{it}) \sim N(\mu_{it}, \tau_{i}) \\ \mu_{it} = \alpha_{i} + \beta_{i} * electrónico_{it} \]
Además: \[ \alpha_{i} \sim N(0, 0.1) \\ \tau_{i} \sim gamma(c_0, c_1) \\ \beta \sim N(0, 0.01) \] donde:
Para comprobar que no existen diferencias significativas entre las mediciones en papel o electrónicas tenemos que: \[ H_{o}: \beta = 0 \\ H_{1}: \beta \neq 0 \]
Siguiendo el modelo utilizando 5000 iteraciones, 3 cadenas y quemando las primeras 1000, se puede observar que cero pertenece al intervalo de la beta, por lo que no se rechaza la hipótesis nula, lo que implica que no existen diferencias entre las mediciones realizadas en papel o electrónicamente.
Además, se puede observar que los valores generales de glucosa como su varianza varían por paciente como se muestra en las gráficas siguientes.
4 horas en las que se miden se refiere a los registros realizados antes del desayuno (08:00), comida (12:00), cena (18:00), y hora de dormir (22:00).outliers. También el primer, segundo y tercer cuartil son más bajos en relación al resto de las mediciones.Con base en lo anterior, sugerimos un modelo jerárquico para responder si existe discrepancia en las mediciones realizadas en las 4 distintas horas del día
El modelo sugerido es
\[ log(y_{it}) \sim N(\mu_{i}, \tau_{i}) \\ \mu_{i} = \alpha_{ij} + \beta_{1} * horario8am_{i} + \beta_{2} * horario12pm_{i} + \beta_{3} * horario18pm_{i} + \beta_{4} * horario22pm_{i} \]
Además:
\[ \alpha_{i} \sim N(0, 0.001) \\ \tau_{i} \sim gamma(c_0, c_1) \\ c0 \sim gamma(1, 1) \\ c1 \sim gamma(1, 1) \\ \beta_i \sim N(0, 0.001) \ ; \ i = 1,2,3,4 \] donde:
El modelo se corrió con 3 distintas cadenas, 5000 iteraciones donde las primeras 1000 son de calentamiento
El resumen estadístico para las betas asociadas a las 4 horas se presentan a continuación y con esto podemos decir que
A continuación se muestra una gráfica boxplot con las betas asociadas a cada momento (desyauno/comida/cena/dormir) y observamos lo que habiamos notado:
Existen individuos muy por debajo de la media global y también por arriba. Incluso, hay algunos que estan por encima del tercer cuartil.
outliers y nuestro rango de datos tiene una gran variabilidad desde 50.04hasta 121.67.Existe discrepancia entre individuos aunque hay algunos
similares entre ellos.
Sí, son más bajas:
33,34,35,36 corresponden a administración de insulina medicada, por lo que ayuda a disminuir el nivel de glucosa en la sangre, por dicha razón estos códigos se asocian con niveles bajos de glucosa.Para esta pregunta se realizó un modelo jeráquico, la jerarquía se aplica a nivel individuo en donde los individuos comparten la precisión. Esto nos permite mejorar la precisión de las estimaciones de los parámetros para los que no se tienen suficientes observaciones. Por ejemplo, en este caso tenemos individuos con menos días observados.
En este primer modelo no se busca predecir el nivel de glucosa sino verificar si existe una diferencia clara entre los días que se ingieren más alimentos de lo usual.
Se crea una base con el promedio diario de glucosa con un indicador que especifica si en ese día el paciente ingirió más alimentos de lo normal.
En esta gráfica se observa que hay una discrepancia clara entre los días que ingieren más alimentos de lo usal. \(Beta1\) es negativa y representa a los días que se consumió una cantidad normal de alimentos y \(Beta2\) es positiva lo cual significa que el nivel de glucosa aumenta si se come más de lo normal.
Se grafica el parámetro \(alpha\) para observar si existe alguna diferencia entre individuos. Observamos que se comportan de manera similar. (Los intervalos se empalman)
mean 2.5% 97.5%
beta.adj[1] -0.03371274 -0.04649043 -0.02045976
beta.adj[2] 0.03371274 0.02045975 0.04649043
[1] 41313.21
[1] 38798.04
Queremos predecir la medición promedio de glucosa para un individuo cuando come una cantidad de alimentos normal o cuando se excede.
Para plantear este modelo, construimos una base con el promedio de las mediciones por individuo, separando por el tipo de consumo.
En exploratorio de los datos vimos que la mediana del nivel de gluscosa es mayor cuando se excede de la cantidad de alimentos habitual.
En la gráfica de ajustados vs osbservados se aprecia que los puntos están alrededor de una linea de 45º. Podemos ver que sobreestimamos mediciones pequeñas y subestimas mediciones grandes.
Los datos observados se encuentran dentro de nuestros intervalos de confianza.
Observamos que las betas indican una discrepancia entre la diferencia de la ingesta de alimentos, pero no es tan discrminante como el modelo anterior. Esto se puede deber a que estamos promediando individuos y no individuos-días.
Se grafica el parámetro \(alpha\) para observar si existe alguna diferencia entre individuos. Observamos que se comportan de forma similar. Sin embargo, en este caso hay alphas cuyos intervalos no se interceptan, entonces estos individuos sí tienen comportamiento distinto en sus mediciones de glucosa.
En este modelo se obtiene un mejor DIC y \(R^2\).
mean 2.5% 97.5%
beta.adj[1] -0.03977104 -0.070395726 -0.007816606
beta.adj[2] 0.03977104 0.007816606 0.070395726
[1] 1022.333
[1] 0.8289517
[1] 915.7123
http://www.stats.ox.ac.uk/~nicholls/MScMCMC15/jags_user_manual.pdf
Doing Bayesian Data Analysis, Second Edition, John K. Kruschke, AP. Bayesian Data Analysis, Andrew Gelman, John B. Carlin, Hal S. Stern and David B. Dunson.
A continuación se muestra el código de los modelos asociados a la pregunta 1,2 y 5.
model
{
#Likelihood
for (i in 1:N) {
y[i] ~ dnorm(mu[i], tau[indiv[i]])
mu[i] <- alpha[indiv[i]] + beta[1] * x1[i]
}
c0 ~ dgamma(.001, 1)
c1 ~ dgamma(1, 1)
# Previas individuos
for(i in 1:nindiv){
alpha[i] ~ dnorm(0, 0.1)
tau[i] ~ dgamma(c0, c1)
}
# Previas horas
beta ~ dnorm(0, 0.01)
## predicciones
for(i in 1:N) {
yf1[i] ~ dnorm(mu[i], tau[indiv[i]])
}
}model
{
#Likelihood
for (i in 1:n) {
y[i] ~ dnorm(mu[i], tau[i])
mu[i] <- alpha[indiv[i]] +
beta[1] * horario_8am[i] +
beta[2] * horario_12pm[i] +
beta[3] * horario_18pm[i] +
beta[4] * horario_22pm[i]
}
# Previas precisión mediciones
for(j in 1:n){
tau[j] ~ dgamma(c0, c1)
}
c0 ~ dgamma(1, 1)
c1 ~ dgamma(1, 1)
# Previas individuos
for(i in 1:n_indiv){
alpha[i] ~ dnorm(0, 0.001)
}
# Previas horas
for(k in 1:m){
beta[k] ~ dnorm(0, 0.001)
}
## estimabilidad
for(k in 1:m){
beta.adj[k] <- beta[k] - mean(beta[])
}
## predicciones
for(i in 1:n) {
yf1[i] ~ dnorm(mu[i], tau[i])
}
}model{
# Verosimilitud
for (i in 1:n) {
y[i] ~ dgamma(r[i], lambda[i])
eta[i] <- alpha[indiv[i]] + beta[x[i]]
r[i] <- lambda[i]*mu[i]
## Liga
log(mu[i]) <- eta[i]
}
#Priors
for (i in 1:n_indiv) {
alpha[i] ~ dnorm(0,tau.a)
}
for (j in 1:2) { beta[j] ~ dnorm(0,0.001) }
tau.a ~ dgamma(0.001,0.001)
for(i in 1:n){lambda[i] ~ dunif(0,1)}
#Prediction 1
for (i in 1:n) { yf1[i] ~ dgamma(r[i], lambda[i]) }
#Estimable quentities
for(j in 1:2){
beta.adj[j] <- beta[j] - mean(beta[])
}
}